Critical Scale-invariance in Healthy Human Heart Rate 
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We demonstrate the robust scale-invariance in the probability density function (PDF) of detrended 
healthy human heart rate increments, which is preserved not only in a quiescent condition, but also in 
a dynamic state where the mean level of heart rate is dramatically changing. This scale-independent 
and fractal structure is markedly different from the scale-dependent PDF evolution observed in a 
turbulent-like, cascade heart rate model. These results strongly support the view that healthy 
human heart rate is controlled to converge continually to a critical state. 

PACS numbers: 87.19.Hh, 05.40.-a, 87.80. Vt, 89.75.Da 



Healthy human heart rate belongs to a special class 
of complex signals, showing long-range temporal correla- 
tions 0,0, non-Gaussianity of the increment's PDF 
and multifractal scaling properties [3J, |4| , and has served 
as the 'benchmark' of choice for studies of biological com- 
plexity. Two alternative mechanisms, both characterized 
by these three features, have been proposed for this heart 
rate complexity: 1) a random (multiplicative) cascade 
process based on the resemblance of the behavior of the 
structure function 5] of heart rate increments to that of 
spatial velocity differences in hydrodynamic turbulence 
[(j; and 2) critical state- like dynamics Q based on the 
resemblance of the scale- invariant properties in heart rate 
to those of many systems operating near the critical point 
of their phase space. To date the exact mechanism for the 
complex heart rate dynamics is unknown. However, as it 
reflects the dynamics of the autonomic nervous system's 
control of heart rate 0, and thus provides potential 
predictors for the mortality of cardiac patients 0, |jj ^| , 
elucidating this mechanism is considered important. 

Lin and Hughson [y] recently reported an analogy 
between turbulence and human heart rate dynamics 
by finding similarity of the structure function — directly 
linked with multifractal formalism — of heart rate in- 
crements to that of spatial velocity differences in a ran- 
dom cascade process proposed as a model of hydrody- 
namic turbulence. One of the common features of such 
cascade-type multifractal models is the evolution in the 
shape of the PDF of the increments from Gaussian at 
large scales to stretched exponential at smaller scales ^| . 
In this study, we disprove this cascade-like assertion and 
demonstrate that heart rate signals do not follow the evo- 
lution in the shape of the (increment) PDF characteristic 
for cascade-like processes, but report for the first time a 
robust scale invariance. As long-range correlation, non- 
Gaussianity and multifractality are also typical charac- 
teristics of a system at the critical point 0, and 



fluctuations in a system at a critical point are generally 
associated with the scale invariance and universal behav- 
ior of the scaling function [pol fl6| , we conclude that such 
robust scale invariance in the increment PDF suggests 
the alternative scenario of the near critical state-like op- 
eration for the healthy heart rate dynamics. 

The long-term heart rate data analyzed have been mea- 
sured as sequential heart interbeat intervals b(i), where 
i is the beat number. We investigate the PDF of heart 
rate increments at different time scales (in beat num- 
bers), where the non-stationarity of the data has been 
eliminated by local detrending [l| . We first integrate the 
b(i), B(m) — Ylj=iHj)i an d the resultant B(m) is di- 
vided into sliding segments of size 2n. Then in each seg- 
ment the best g-th order polynomial is fit to the data. 
The differences A s B(i) = B*(i + s) — B*(i) at a scale s 
are obtained by sliding in time over the segments, where 
B*{i) is a deviation from the polynomial fit. By this 
procedure, the (q — l)-th order polynomial trends are 
eliminated and we analyze the whole PDF of A s B(i). 

Using this method, we analyze two experimental and 
two computer-generated data sets. The first data set con- 
sists of daytime (12:00-18:00 hrs) heart rate data with a 
length of up to 4 x 10 4 heartbeats from 50 healthy sub- 
jects (10 females and 40 males; ages 21 — 76 years) with- 
out any known disease affecting the autonomic control of 
heart rate [Fig. da)]. Details of the recruitment of the 
subjects, screening for medical problems, protocols and 
the data collection are described in Sakata et ai. |l7|. 
The data were collected during normal daily life. The 
second experimental data set consists of data of seven 26- 
hour long periods (up to 10 5 beats), collected when the 
subjects (7 males; ages 21—30 years) underwent 'constant 
routine' (CR) protocol, where known behavioral factors 
affecting heart rate (e.g. exercise, diet, postural changes 
and sleep) are eliminated [Fig. IHb)] j4Lll8l|. 

In order to test the possible presence of non-linear 
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FIG. 1: (a) A representative record of daytime (12:00-18:00 
hrs) heart interbeat intervals for a healthy subject, (b) An 
example of heartbeat intervals for a healthy subject during 
constant routine protocol P. Il8j|. (c) The surrogate data for 
(b). (d) Data generated by a cascade heart rate increment 
model [g. The parameters used for simulation are J — 15, 
Rt = 2 and = 7? t " 2 - 5 ~ j/J (see Ref. @] for details). 

mechanisms in complex heart rate dynamics, we a pply 
the surrogate data test to the CR protocol data |19| . 
We generate a surrogate data set with the same Fourier 
amplitudes and distributions as the original increments 
in the CR protocol data [Fig. (He)]. Since only linear 
temporal correlation of b(i + 1) — b(i) is retained in the 
surrogate data, a comparison with the raw data can be 
used to test whether the PDF of 'velocity' increments 
A s B(i) possesses some non-linear mechanism inherent to 
it. Finally, we generate heart rate increments of compara- 
ble lengths, following the 'turbulence-like' scenario from 
the random cascade model proposed recently by Lin and 
Hughson [Fig. Old)] @. 

PDF's of A s B(i) for healthy humans, which are stan- 
dardized by dividing the heart rate increments in each 
record by the standard deviation, are non-Gaussian in 
shape for a wide range of scales 8 < s < 4, 096 irrespec- 
tive of whether the subjects were in their normal daily 



routine [Fig. 0[a)] or in CR [Fig. Efb)]. In contrast, the 
PDF's of the surrogate data are near Gaussian, although 
non-Gaussianity with the fat tails close to those in the 
observed data is still encountered at fine scales [Fig. Etc)]. 
The difference between the healthy human and surrogate 
data indicates that the observed non-Gaussian behavior 
is related to non-linear features of the healthy heart rate 
dynamics. The PDF's of the cascade model show contin- 
uous deformation and the appearance of fat tails when 
going from large to small scales [Fig. Eld)] ■ 

For a quantitative comparison, we fit the data to the 
following function based on Castaing's equation |l2^ : 

P.(x) = J P L (^) ^G SiL Qjia)d(lna), 

where Pl is the increment PDF at a large scale L > s, 
and the self-similarity kernel G s ,l determines the nature 
of the cascade- type multiplicative process. Here we as- 
sume Pl and G s ,l are both Gaussian, 

and investigate the scale dependence of A 2 . The fit of the 
PDF of actual heart rate increments to Castaing's equa- 
tion is indeed almost perfect, especially within ±3 times 
standard deviation, even for a single record [Fig. EJa)], 
and robust in terms of the effect of the order of defend- 
ing polynomials on the estimation of A 2 , if the order is 
greater than two [Fig. EJb)] . In the following, we use the 
third order detrending for the estimation of A 2 . 

Within the turbulent cascade picture, the parameter 
A 2 can be interpreted as being proportional to the num- 
ber of cascade steps and is known to decrease linearly 
with log s 0, |2(j ■ The cascade heart rate model studied 
here also shows this effect [Fig.[3fc)]. In contrast, the 
scale dependency of A 2 for healthy heart rate increments 
is remarkably different [Fig.^c)]. Especially during CR, 
we cannot see any decrease in A 2 with logs. There is 
no significant difference in the average A 2 at different 
scales, tested by the analysis of variance, over the range 
of 23-2,048 beats for healthy subjects during daily rou- 
tine and 8-4,096 beats during CR [F(13, 686) = 1.73 and 
F(18, 114) = 1.70, respectively, p > 0.05], and the slopes 
of A 2 vs. log s are much closer to zero, which means the 
absence of cascade steps across the scales in the corre- 
sponding range. 

In addition, when the PDF's at different scales are 
superimposed [Fig. EJd)], all the data collapse on the 
same curve, which is one of the characteristic features ob- 
served in fluctuations at a critical point |2l|. The range 
of scales where this scale-invariance of PDF is observed, 
spanning from about 10 beats to a few thousand heart- 
beats, is compatible with that of the robust, behavioral- 
independent 1// scaling 0] and multifractality 0] of 
heart rate. The scale-invariance in the PDF is also ro- 
bust in the sense that it is observed not only during CR 
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FIG. 2: Deformation of increment PDF's across scales. Standardized PDF's (in logarithmic scale) of A s B(i) for different time 
scales are shown for (from top to bottom) s = 8, 16, 32, 64, 128, 256, 512, 1024 beats. These PDF's are estimated from all 
samples in each group. The dashed line is a Gaussian PDF for comparison, (a) The PDF's from daytime (12:00-18:00 hrs) 
heart rate time series from healthy subjects, (b) From healthy subjects during constant routine protocol, (c) From surrogate 
time series for (b). (d) From a cascade model. Note in the last case (d) the continuous deformation and the appearance of fat 
tails when going from large to fine scales. In solid lines, we superimpose the deformation of the PDF using Castaing's equation 
with the log-normal self-similarity kernel, providing an excellent fit to the data. 



but also during normal daily life, where behavioral mod- 
ifiers of heart rate dramatically change the mean level of 
heart rate [e.g. Fig.^a)]. We thus find a novel property 
of scale invariance in healthy human heart rate dynamics, 
reminiscent of systems in a critical state. In particular, 
the invariance discovered strongly supports the view that 
the healthy human heart is controlled to converge contin- 
ually to a critical state. Such a critical point itself may, 
however, be shifted by the effects of the external and/or 
internal environment. 

Struzik et al. recently demonstrated that modify- 
ing the relative importance of either the sympathetic or 
the parasympathetic branch of the autonomic nervous 
system leads to a substantial decrease in 1// scaling, 
showing that 1// scaling in healthy heart rate requires 
the existence of and the intricate balance between the an- 
tagonistic activity of these two branches. They further 
suggest the view of cardiac neuroregulation as a system 
in a critical state |23j . and permanently out of equilib- 
rium, in which concerted interplay of the sympathetic 
and parasympathetic nervous systems is required for pre- 
serving momentary 'balance'. Our findings provide more 
direct evidence for this. The precise mechanism respon- 
sible for critical heart rate dynamics requires further re- 
search. It is of note, however, that there exists a physi- 
ological model for the dynamics of cardiac neuroregula- 



tion |2J], equipped with antagonistic and multiplicative 
delayed feedback loops, within time scales where the crit- 
ical scale-invariance in heart rate is observed. The mech- 
anism of the critical mode of operation could be clarified 
by investigating essential dynamics in such a 'first prin- 
ciples', non-linear physiological model. 

The functional advantage of the heart rate control 
system being in a critical state remains an open ques- 
tion. However, an analogy with other critical phenom- 
ena might help to understand this. Studies on trans- 
port properties through complex networks |2ol l2o| have 
demonstrated maximum efficiency of transportation at 
the critical point, which is the phase transition point from 
an 'uncrowded' state to a 'congested' state in the trans- 
portation routes. Thus, our results may indicate that 
the central neuroregulation continually brings the heart 
to a critical state to maximize its functional ability, cop- 
ing with the continually changing pre- and after-load on 
the heart. This may be particularly important in under- 
standing the widely reported evidence that decreased 1/f 
variability 0, [|| , especially in the low frequency region 
0, 0| , is associated with increased mortality in cardiac 
patients. To date there have been no successful attempts 
to provide a satisfactory explanation for this. We sug- 
gest for the first time that a breakdown of the optimal 
control achieved in the critical state might be related to 
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FIG. 3: (a) Standardized PDF for a single subject during constant routine protocol [see Fig. 0b)]. In solid lines, we 
superimpose the PDF's using Castaing's equation, (b) Dependence of the fitting parameter A 2 of Castaing's equation on the 
order of detrending polynomials, (c) Dependence of the fitting parameter A 2 on the scale s. The error bars indicate the standard 
error of the group averages, (d) Superposition of standardized PDF's at different scales shown in Fig. and|!J>, where we 
use the scale range 8 < s < 2048 and 8 < s < 4096, respectively. In solid lines, we superimpose the PDF using the Castaing's 
equation (A 2 = 0.16 for both groups). 



this clinically important phenomenon. 
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